FC Development Main Figures

schaefer200.parcel.labels <- read.csv("/cbica/projects/network_replication/atlases/parcellations/schaefer200x7_regionlist_final.csv")

schaefer200_SAaxis <- read.csv("/cbica/projects/network_replication/SAaxis/schaefer200x7_SAaxis.csv")
schaefer200_SAaxis <- schaefer200_SAaxis %>% rename(region=label)

Figure 1

# Function for combining final GAM output dfs for figures
# @param metric A character string of connectivity metric (i.e. "GBC")
# @param atlas A character string of atlas name
combineGAMdfs <- function(metric, dataset){
  # Combine into Final Dfs
  gam_output <- get(paste0("gam.", metric, ".age.", dataset))
  SAaxis <- schaefer200_SAaxis
  df.list <- list(SAaxis, gam_output)
  axis <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
  axis$SA.axis_rank <- as.numeric(axis$SA.axis_rank)
  return(axis)
}

# Function for calculating number of significant parcels (FDR corrected)
# @param axis A dataframe with SA-axis rank and GAM results
sig_parcels <- function(axis) {
  axis$Anova.age.pvalue.fdr <- p.adjust(axis$Anova.age.pvalue, method=c("fdr"))
  #cat(sprintf("There are %s/%s significant parcels", sum(axis$Anova.age.pvalue.fdr < 0.05, na.rm=TRUE), nrow(axis)))
  
  axis$significant.fdr <- axis$Anova.age.pvalue.fdr < 0.05
  axis$significant.fdr[axis$significant.fdr == TRUE] <- 1
  axis$significant.fdr[axis$significant.fdr == FALSE] <- 0
  return(axis)
}

gam.GBC.age.PNC <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMresults.%3$s.age.%4$s.csv", derivs_root, "PNC", "GBC", "schaefer200"))
gam.GBC.age.PNC$label <- gsub("Network", "7Network", gam.GBC.age.PNC$label)
gam.GBC.age.PNC <- gam.GBC.age.PNC %>% rename(region=label)

gam.GBC.age.NKI <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMresults.%3$s.age.%4$s.csv", derivs_root, "NKI", "GBC", "schaefer200"))
gam.GBC.age.NKI$label <- gsub("Network", "7Network", gam.GBC.age.NKI$label)
gam.GBC.age.NKI <- gam.GBC.age.NKI %>% rename(region=label)

gam.GBC.age.HCPD <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMresults.%3$s.age.%4$s.csv", derivs_root, "HCPD", "GBC", "schaefer200"))
gam.GBC.age.HCPD$label <- gsub("Network", "7Network", gam.GBC.age.HCPD$label)
gam.GBC.age.HCPD <- gam.GBC.age.HCPD %>% rename(region=label)

gam.GBC.age.HBN <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMresults.%3$s.age.%4$s.csv", derivs_root, "HBN", "GBC", "schaefer200"))
gam.GBC.age.HBN$label <- gsub("Network", "7Network", gam.GBC.age.HBN$label)
gam.GBC.age.HBN <- gam.GBC.age.HBN %>% rename(region=label)
 
GBC.axis.PNC <- combineGAMdfs("GBC", "PNC") %>% sig_parcels
GBC.axis.NKI <- combineGAMdfs("GBC", "NKI") %>% sig_parcels
GBC.axis.HCPD <- combineGAMdfs("GBC", "HCPD") %>% sig_parcels
GBC.axis.HBN <- combineGAMdfs("GBC", "HBN") %>% sig_parcels

Figure 1a

# spin test
perm.id.full_schaefer200 <- readRDS("/cbica/projects/network_replication/software/rotate_parcellation/schaefer200x7.coords_sphericalrotations_N10k_seed10.rds")
 
# GBC - make df.dev.spin_GBC_schaefer200 
parcel.labels <- schaefer200.parcel.labels$label
df.dev.spin_GBC.axis.PNC <- rbind(GBC.axis.PNC[1:c(length(parcel.labels)/2),], GBC.axis.PNC[c(length(parcel.labels)/2+1):length(parcel.labels),]) 
df.dev.spin_GBC.axis.NKI <- rbind(GBC.axis.NKI[1:c(length(parcel.labels)/2),], GBC.axis.NKI[c(length(parcel.labels)/2+1):length(parcel.labels),])
df.dev.spin_GBC.axis.HCPD <- rbind(GBC.axis.HCPD[1:c(length(parcel.labels)/2),], GBC.axis.HCPD[c(length(parcel.labels)/2+1):length(parcel.labels),])
df.dev.spin_GBC.axis.HBN <- rbind(GBC.axis.HBN[1:c(length(parcel.labels)/2),], GBC.axis.HBN[c(length(parcel.labels)/2+1):length(parcel.labels),])

# spatial correlation of GAM adjusted age Rsq between all pairs of datasets
corr.GBC.axis.PNC_GBC.axis.NKI <- cor.test(GBC.axis.PNC$GAM.age.AdjRsq, GBC.axis.NKI$GAM.age.AdjRsq, method=c("pearson"), exact=F)$estimate # r = 0.8804216 
corr.GBC.axis.PNC_GBC.axis.HCPD <- cor.test(GBC.axis.PNC$GAM.age.AdjRsq, GBC.axis.HCPD$GAM.age.AdjRsq, method=c("pearson"), exact=F)$estimate # r = 0.8279745
corr.GBC.axis.PNC_GBC.axis.HBN <- cor.test(GBC.axis.PNC$GAM.age.AdjRsq, GBC.axis.HBN$GAM.age.AdjRsq, method=c("pearson"), exact=F)$estimate # r = 0.629025  
corr.GBC.axis.NKI_GBC.axis.HCPD <- cor.test(GBC.axis.HCPD$GAM.age.AdjRsq, GBC.axis.NKI$GAM.age.AdjRsq, method=c("pearson"), exact=F)$estimate # r = 0.6926665 
corr.GBC.axis.NKI_GBC.axis.HBN <- cor.test(GBC.axis.HBN$GAM.age.AdjRsq, GBC.axis.NKI$GAM.age.AdjRsq, method=c("pearson"), exact=F)$estimate # r = 0.4665967
corr.GBC.axis.HBN_GBC.axis.HCPD <- cor.test(GBC.axis.HBN$GAM.age.AdjRsq, GBC.axis.HCPD$GAM.age.AdjRsq, method=c("pearson"), exact=F)$estimate # r = 0.7534226 
 
# do spin test for each pair of datasets' spatial correlation 
p_spin_GBC.axis.PNC_GBC.axis.NKI  <- perm.sphere.p(GBC.axis.PNC$GAM.age.AdjRsq, GBC.axis.NKI$GAM.age.AdjRsq, perm.id.full_schaefer200, corr.type='pearson')  
p_spin_GBC.axis.PNC_GBC.axis.HCPD <- perm.sphere.p(GBC.axis.PNC$GAM.age.AdjRsq, GBC.axis.HCPD$GAM.age.AdjRsq, perm.id.full_schaefer200, corr.type='pearson')  
p_spin_GBC.axis.PNC_GBC.axis.HBN  <- perm.sphere.p(GBC.axis.PNC$GAM.age.AdjRsq, GBC.axis.HBN$GAM.age.AdjRsq, perm.id.full_schaefer200, corr.type='pearson')  
p_spin_GBC.axis.NKI_GBC.axis.HCPD  <- perm.sphere.p(GBC.axis.HCPD$GAM.age.AdjRsq, GBC.axis.NKI$GAM.age.AdjRsq, perm.id.full_schaefer200, corr.type='pearson')  
p_spin_GBC.axis.NKI_GBC.axis.HBN  <- perm.sphere.p(GBC.axis.HBN$GAM.age.AdjRsq, GBC.axis.NKI$GAM.age.AdjRsq, perm.id.full_schaefer200, corr.type='pearson')  
p_spin_GBC.axis.HBN_GBC.axis.HCPD  <- perm.sphere.p(GBC.axis.HBN$GAM.age.AdjRsq, GBC.axis.HCPD$GAM.age.AdjRsq, perm.id.full_schaefer200, corr.type='pearson')  

# make correlation matrix
dat2 <- bind_cols(GBC.axis.PNC$region, as.numeric(GBC.axis.PNC$GAM.age.AdjRsq), GBC.axis.NKI$GAM.age.AdjRsq, GBC.axis.HCPD$GAM.age.AdjRsq, GBC.axis.HBN$GAM.age.AdjRsq)
names(dat2) <- c("label", "PNC", "NKI", "HCP-D", "HBN")
corr <- round(cor(dat2[,-1], use = "complete.obs"), 5)
  
# compute matrix of correlation p-values
p.mat <- cor_pmat(dat2[,-1])

 
# Visualize the correlation matrix
correlation_matrix <- ggcorrplot(corr,  
   lab = TRUE, lab_col = "white", type = "lower",
   outline.col = "white",
   ggtheme = ggplot2::theme_void, lab_size = 8, 
   hc.order = FALSE, insig = "blank", p.mat=p.mat)  +
   scale_fill_gradient2(low = "white",
  mid = "#E5B9ADFF",
  high = "#A50026FF", midpoint=0.5,limit= c(0, 1)) + 
 theme(legend.text = element_text(size = 18), 
       legend.title = element_blank(), 
       legend.key.height = unit(1, 'cm'),
       legend.key.width = unit(2, 'cm'),
       legend.position = "bottom",
       legend.margin=margin(10,0,0,0),
       axis.text=element_text(size=22, color = "black"),
       axis.text.y = element_text(hjust = 0, size = 22),
       axis.text.x = element_text(hjust = 0.5, angle=0, size = 22),
       plot.margin = margin(1, 1, 1, 1, "cm"))


 

# plot figure 1a
grid.arrange(arrangeGrob(correlation_matrix, bottom = textGrob("Correlation", 
                                                     gp=gpar(col="black", fontsize=22), x = 0.6, y = 1.5)))

write.csv(corr, "/cbica/projects/network_replication/manuscript/figure_data/figure1a_corrplot.csv")

Figure 1b

figure1b <- ggplot() + geom_brain(data=schaefer200_SAaxis, 
                        atlas=schaefer7_200, 
                        mapping=aes(fill=SA.axis_rank),
                        show.legend=TRUE, 
                        position = position_brain(c('right lateral', 'right medial'))) +
  paletteer::scale_fill_paletteer_c("grDevices::RdYlBu", direction=-1) +
  theme_void() +
  theme(legend.text = element_text(size = 18), 
       legend.title = element_blank(), 
       legend.key.height = unit(1, 'cm'),
       legend.key.width = unit(2.3, 'cm'),
       legend.position = "bottom",
       legend.margin=margin(10,0,0,0),
       plot.margin = margin(0, 0, 0, 0, "cm"))
  


# plot figure 1b
grid.arrange(arrangeGrob(figure1b, bottom = textGrob("Sensorimotor-Association Axis Rank", 
                                                     gp=gpar(col="black", fontsize=20), x = 0.5, y = 4.6)))

write.csv(schaefer200_SAaxis, "/cbica/projects/network_replication/manuscript/figure_data/figure1b_SAaxis.csv")

Figure 1c

PNC <- GBC.axis.PNC
NKI <- GBC.axis.NKI
HCPD <- GBC.axis.HCPD
HBN <- GBC.axis.HBN

PNC <- PNC %>% mutate(SA.axis_rank_signif = ifelse(significant.fdr == 1, SA.axis_rank, NA))
PNC <- bind_cols(PNC$region, PNC$GAM.age.AdjRsq, PNC$significant.fdr)
names(PNC) <- c("region", "GAM.age.AdjRsq", "significant.fdr")
PNC$significant.fdr <- gsub("0", "Non_Sig", PNC$significant.fdr)
PNC$significant.fdr <- gsub("1", "Sig", PNC$significant.fdr)
 
  
 
NKI <- NKI %>% mutate(SA.axis_rank_signif = ifelse(significant.fdr == 1, SA.axis_rank, NA))
NKI <- bind_cols(NKI$region, NKI$GAM.age.AdjRsq, NKI$significant.fdr)
names(NKI) <- c("region", "GAM.age.AdjRsq", "significant.fdr")
NKI$significant.fdr <- gsub("0", "Non_Sig", NKI$significant.fdr)
NKI$significant.fdr <- gsub("1", "Sig", NKI$significant.fdr)

 


HCPD <- HCPD %>% mutate(SA.axis_rank_signif = ifelse(significant.fdr == 1, SA.axis_rank, NA))
HCPD <- bind_cols(HCPD$region, HCPD$GAM.age.AdjRsq, HCPD$significant.fdr)
names(HCPD) <- c("region", "GAM.age.AdjRsq", "significant.fdr")
HCPD$significant.fdr <- gsub("0", "Non_Sig", HCPD$significant.fdr)
HCPD$significant.fdr <- gsub("1", "Sig", HCPD$significant.fdr)

HBN <- HBN %>% mutate(SA.axis_rank_signif = ifelse(significant.fdr == 1, SA.axis_rank, NA))
HBN <- bind_cols(HBN$region, HBN$GAM.age.AdjRsq, HBN$significant.fdr)
names(HBN) <- c("region", "GAM.age.AdjRsq", "significant.fdr")
HBN$significant.fdr <- gsub("0", "Non_Sig", HBN$significant.fdr)
HBN$significant.fdr <- gsub("1", "Sig", HBN$significant.fdr)
 
 

ageEffect_surf.fig <- function(dataset_name, title){
  surf.fig <- ggplot() + geom_brain(data=get(dataset_name), 
                        atlas=schaefer7_200, 
                        mapping=aes(fill=GAM.age.AdjRsq, colour=significant.fdr, size=I(0.3)),
                        show.legend=TRUE, 
                        position = position_brain(hemi~side)) + scale_colour_manual(values = c(Sig = "black", Non_Sig = "grey84")) + 
      scale_fill_gradientn(colors= c("#8a0f63", "#FFFFFF","#f4b100"), limits = c(-0.05,0.15),  oob=squish,
                           values=rescale(c(-0.05,0,0.15))) +  theme_void() + 
    labs(title=title) + 
    guides(color=FALSE) + 
    theme(legend.position = "bottom",
        legend.key.height = unit(1, 'cm'),
        legend.key.width = unit(2.3, 'cm'),
        legend.margin=margin(10,0,0,0),
        legend.text = element_text(size=24),
        legend.title = element_blank(),
        plot.title = element_text(size=20, hjust = 0.5)) 
     
  return(surf.fig)
}
  
 

ageEffect_surf_PNC <- ageEffect_surf.fig("PNC", "PNC (N = 1207)")
ageEffect_surf_NKI <- ageEffect_surf.fig("NKI", "NKI (N = 397)")
ageEffect_surf_HCPD <- ageEffect_surf.fig("HCPD", "HCP-D (N = 625)")
ageEffect_surf_HBN <- ageEffect_surf.fig("HBN", "HBN (N = 1126)")

figure1c_alldatasets <- ggarrange(ageEffect_surf_PNC, ageEffect_surf_HCPD, ageEffect_surf_NKI, ageEffect_surf_HBN, common.legend=TRUE, legend = c('bottom') )

x.grob <- textGrob(expression(paste("Age Effect (Delta Adj", " R"^2, ")")), 
                   gp=gpar(col="black", fontsize=20))
 

grid.arrange(arrangeGrob(figure1c_alldatasets, bottom = x.grob))

write.csv(PNC, "/cbica/projects/network_replication/manuscript/figure_data/figure1c_PNC_age_effect.csv")
write.csv(NKI, "/cbica/projects/network_replication/manuscript/figure_data/figure1c_NKI_age_effect.csv")
write.csv(HCPD, "/cbica/projects/network_replication/manuscript/figure_data/figure1c_HCPD_age_effect.csv")
write.csv(HBN, "/cbica/projects/network_replication/manuscript/figure_data/figure1c_HBN_age_effect.csv")

Figure 2

## FIGURE: developmental trajectories, centered on zero
# @param atlas A character string of atlas name 
# @param metric A character string of connectivity metric (i.e. "GBC")
make_smooths_fig_centered <- function(dataset){
  smooth_fits <- get(paste0(dataset, "_smooths"))
  smooths_fig_centered <- ggplot(smooth_fits,aes(age,est,group=SA.axis_rank)) + 
    geom_line(data = smooth_fits, size=.7, alpha = .6, aes(color=SA.axis_rank)) + 
    ylim(-0.042, 0.053) + xlim(5, 23) +
    paletteer::scale_color_paletteer_c("grDevices::RdYlBu", direction = -1, limits = c(min(smooth_fits$SA.axis_rank), max(smooth_fits$SA.axis_rank)), oob = squish) + 
    theme(
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      axis.line = element_line(color = "black"),
      axis.text=element_text(size=28, color = "black"),
      panel.background=element_blank(),  
      legend.position = "none", 
      plot.margin = margin(t = 1, r = 0, b = , l = 1, unit = "cm")) 
  return(smooths_fig_centered)
}

# subsets gam.smooths (long_df) to extract specific parcels and creates a new dataframe (to make dev trajectory figures)
# @param SA_parcelnames, A character string indicating which parcels to plot (i.e "SM", "default", etc.; may be different for different atlases)
smooths_repParcels <- function(SA_ranks, dataset){
  gam.smooths.centered <- get(paste0(dataset, "_smooths"))
  gam.smooths.ordered <- gam.smooths.centered[order(gam.smooths.centered$SA.axis_rank),]
  parcel_rows <- which(gam.smooths.ordered$SA.axis_rank %in% SA_ranks) # trying to extract all rows that match SA_ranks
  repParcels <- gam.smooths.ordered[parcel_rows,]
  return(repParcels)
}

## FIGURE: developmental trajectories for representative parcels
# @param dataset A character string for dataset
# @param parcels A character string of parcel(s) whose trajectory you want to plot (i.e. "SomMot", "visual")
# @param color_hexcode A character string for color of geom_line
make_repParcel.fig <- function(dataset, parcels, color_hexcode) {
  df <- get(paste0("repParcels_", dataset))
  smooth_fits <- df[[parcels]]
  main.plot <- ggplot(smooth_fits,aes(age,est,group=SA.axis_rank)) + 
    geom_line(data = smooth_fits, size=1.5, alpha = .6,  colour=color_hexcode) + ylim(-0.042, 0.053) + xlim(5, 23) +
    theme(
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      axis.line = element_line(color = "black"),
      axis.text=element_text(size=28, color = "black"),
      panel.background=element_blank(),  
      legend.position = "none", plot.margin = margin(t = 1, r = 0, b = , l = 1, unit = "cm"))  
  return(main.plot)
}

PNC_smooths <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMsmoothfits.%3$s.age.%4$s_centered.csv", derivs_root, "PNC", "GBC", "schaefer200"))
NKI_smooths <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMsmoothfits.%3$s.age.%4$s_centered.csv", derivs_root, "NKI", "GBC", "schaefer200"))
HCPD_smooths <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMsmoothfits.%3$s.age.%4$s_centered.csv", derivs_root, "HCPD", "GBC", "schaefer200"))
HBN_smooths <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/GAMsmoothfits.%3$s.age.%4$s_centered.csv", derivs_root, "HBN", "GBC", "schaefer200"))

devTraj.centered.PNC <- make_smooths_fig_centered("PNC")
devTraj.centered.NKI <- make_smooths_fig_centered("NKI")
devTraj.centered.HCPD <- make_smooths_fig_centered("HCPD")
devTraj.centered.HBN <- make_smooths_fig_centered("HBN")

SomMot <- which(str_detect(schaefer200_SAaxis$region, "SomMot")) # mean SA rank = 40
SomMot_ranks <- schaefer200_SAaxis$SA.axis_rank[SomMot]  

attention <- which(str_detect(schaefer200_SAaxis$region, "SalVentAttn")) # mean SA rank = 113
attention_ranks <- schaefer200_SAaxis$SA.axis_rank[attention]  
 
default <- which(str_detect(schaefer200_SAaxis$region, "Default")) # mean SA rank = 155
default_ranks <- schaefer200_SAaxis$SA.axis_rank[default]  
 

schaefer200_parcelRanks <- list(SomMot_ranks, attention_ranks, default_ranks)
repParcels_PNC <- lapply(schaefer200_parcelRanks, smooths_repParcels, dataset="PNC") 
names(repParcels_PNC) <- c("SomMot", "attention", "default")

repParcels_NKI <- lapply(schaefer200_parcelRanks, smooths_repParcels, dataset="NKI")  
names(repParcels_NKI) <- c("SomMot", "attention", "default")

repParcels_HCPD <- lapply(schaefer200_parcelRanks, smooths_repParcels, dataset="HCPD")  
names(repParcels_HCPD) <- c("SomMot", "attention", "default")

repParcels_HBN <- lapply(schaefer200_parcelRanks, smooths_repParcels, dataset="HBN")  
names(repParcels_HBN) <- c("SomMot", "attention", "default")
  
## set up dataframe for making figures
SAaxis_repParcels_schaefer200 <- schaefer200_SAaxis
SAaxis_repParcels_schaefer200$SA.axis_rank <- as.numeric(SAaxis_repParcels_schaefer200$SA.axis_rank)
names(SAaxis_repParcels_schaefer200)[3] <- "region"
schaefer200_SAaxis_fig.df <- SAaxis_repParcels_schaefer200 %>% mutate(SomMot = ifelse(SA.axis_rank %in% SomMot_ranks, 1, 0)) %>% 
  mutate(default = ifelse(SA.axis_rank %in% default_ranks, 1, 0))  %>%  
  mutate(attention = ifelse(SA.axis_rank %in% attention_ranks, 1, 0)) 

PNC_sensorimotor_figs <- make_repParcel.fig("SomMot", dataset="PNC", color_hexcode="#24A6A8FF")
PNC_midaxis_figs <- make_repParcel.fig("attention", dataset="PNC", color_hexcode="#EEAC32FF")
PNC_assoc_figs <- make_repParcel.fig("default", dataset="PNC", color_hexcode="#C73000FF")

NKI_sensorimotor_figs <- make_repParcel.fig("SomMot", dataset="NKI", color_hexcode="#24A6A8FF")
NKI_midaxis_figs <- make_repParcel.fig("attention", dataset="NKI", color_hexcode="#EEAC32FF")
NKI_assoc_figs <- make_repParcel.fig("default", dataset="NKI", color_hexcode="#C73000FF")

HCPD_sensorimotor_figs <- make_repParcel.fig("SomMot", dataset="HCPD", color_hexcode="#24A6A8FF")
HCPD_midaxis_figs <- make_repParcel.fig("attention", dataset="HCPD", color_hexcode="#EEAC32FF")
HCPD_assoc_figs <- make_repParcel.fig("default", dataset="HCPD", color_hexcode="#C73000FF")

HBN_sensorimotor_figs <- make_repParcel.fig("SomMot", dataset="HBN", color_hexcode="#24A6A8FF")
HBN_midaxis_figs <- make_repParcel.fig("attention", dataset="HBN", color_hexcode="#EEAC32FF")
HBN_assoc_figs <- make_repParcel.fig("default", dataset="HBN", color_hexcode="#C73000FF")

figure2_alldatasets <- plot_grid(devTraj.centered.PNC, 
                               PNC_sensorimotor_figs + labs(title = "Somatomotor") + 
                                 theme(plot.title=element_text(size=28, face="bold", hjust=0.5, vjust=5)), 
                               PNC_midaxis_figs + labs(title = "Attention") + 
                                 theme(plot.title=element_text(size=28, face="bold", hjust=0.5, vjust=5)), 
                               PNC_assoc_figs + labs(title = "Default") + 
                                 theme(plot.title=element_text(size=28, face="bold", hjust=0.5, vjust=5)), 
          devTraj.centered.NKI, NKI_sensorimotor_figs,  NKI_midaxis_figs, NKI_assoc_figs,
          devTraj.centered.HCPD, HCPD_sensorimotor_figs, HCPD_midaxis_figs, HCPD_assoc_figs, 
          devTraj.centered.HBN, HBN_sensorimotor_figs, HBN_midaxis_figs, HBN_assoc_figs, ncol=4, 
          labels=c("a", "e", "i", "m", "b", "f", "j", "n","c", "g", "k", "o", "d", "h", "l", "p"), 
          label_size=28, label_fontface = 'bold', label_y=rep(1,32), label_x = 0.01, byrow=TRUE)



x.grob <- textGrob("Age (years)", 
                   gp=gpar(fontface="bold", col="black", fontsize=28))
y.grob <- textGrob("Functional Connectivity Strength (Zero-Centered)", 
                   gp=gpar(fontface="bold", col="black", fontsize=28), rot=90)

grid.arrange(arrangeGrob(figure2_alldatasets, left = y.grob, bottom = x.grob))

# row 1 = PNC
# row 2 = NKI
# row 3 = HCP-D
# row 4 = HBN


write.csv(PNC_smooths, "/cbica/projects/network_replication/manuscript/figure_data/figure2aeim_PNC_devtrajectories.csv")
write.csv(NKI_smooths, "/cbica/projects/network_replication/manuscript/figure_data/figure2bfjn_NKI_devtrajectories.csv")
write.csv(HCPD_smooths, "/cbica/projects/network_replication/manuscript/figure_data/figure2cgko_HCPD_devtrajectories.csv")
write.csv(HBN_smooths, "/cbica/projects/network_replication/manuscript/figure_data/figure2dhlp_HBN_devtrajectories.csv")

Figure 3

# function for creating the final df.dev.spin dataframe for permutation testing
prepDfSpin <- function(atlas, dataset, metric) {
  parcel.labels <- get(paste0(atlas, ".parcel.labels"))
  parcel.labels <- parcel.labels$label
  axis.df <- get(paste0(metric, ".axis.", dataset))
  figureDF <- axis.df
  if(atlas=="gordon"){
    df.dev.spin <- rbind(figureDF[1:161,], figureDF[162:333,]) #format df as left hemisphere -> right hemisphere for spin tests
  } else if(str_detect(atlas, "schaefer") | str_detect(atlas, "glasser")) {
    df.dev.spin <- rbind(figureDF[1:c(length(parcel.labels)/2),], figureDF[c(length(parcel.labels)/2+1):length(parcel.labels),]) #format df as left hemisphere -> right hemisphere for spin tests
  }  
  return(df.dev.spin)
}

# make df.dev.spin.<dataset> (formatted df's for applying spin test)
df.dev.spin.PNC <- prepDfSpin("schaefer200","PNC", "GBC")
df.dev.spin.NKI <- prepDfSpin("schaefer200","NKI", "GBC")
df.dev.spin.HCPD <- prepDfSpin("schaefer200","HCPD", "GBC")
df.dev.spin.HBN <- prepDfSpin("schaefer200","HBN", "GBC")

# compute spatial correlation between age effect and S-A axis rank
PNC_corr <- round(cor.test(GBC.axis.PNC$GAM.age.AdjRsq, as.numeric(GBC.axis.PNC$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
NKI_corr <- round(cor.test(GBC.axis.NKI$GAM.age.AdjRsq, as.numeric(GBC.axis.NKI$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
HCPD_corr <- round(cor.test(GBC.axis.HCPD$GAM.age.AdjRsq, as.numeric(GBC.axis.HCPD$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
HBN_corr <- round(cor.test(GBC.axis.HBN$GAM.age.AdjRsq, as.numeric(GBC.axis.HBN$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)

#Load Spin Test Parcel Rotation Matrix
perm.id.full_schaefer200 <- readRDS("/cbica/projects/network_replication/software/rotate_parcellation/schaefer200x7.coords_sphericalrotations_N10k.rds")



# spin permutation tests:
PNC_pspin <- perm.sphere.p(as.numeric(df.dev.spin.PNC$GAM.age.AdjRsq), as.numeric(df.dev.spin.PNC$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

NKI_pspin <- perm.sphere.p(as.numeric(df.dev.spin.NKI$GAM.age.AdjRsq), as.numeric(df.dev.spin.NKI$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

HCPD_pspin <- perm.sphere.p(as.numeric(df.dev.spin.HCPD$GAM.age.AdjRsq), as.numeric(df.dev.spin.HCPD$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman') 

HBN_pspin <- perm.sphere.p(as.numeric(df.dev.spin.HBN$GAM.age.AdjRsq), as.numeric(df.dev.spin.HBN$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman') 

corr_pspin <- data.frame(cbind(c(PNC_corr, NKI_corr, HCPD_corr, HBN_corr), c(PNC_pspin, NKI_pspin, HCPD_pspin, HBN_pspin))) %>% setNames(c("rho", "pspin"))
 
rownames(corr_pspin) <- c("PNC", "NKI", "HCPD", "HBN")

gam_figure3_outliers <- function(dataset, axis, metric, atlas, r_text, x_pos, y_pos, ylim_lower, ylim_upper){
  axis <- axis %>% mutate(SA.axis_rank_signif = ifelse(significant.fdr == 1, SA.axis_rank, NA))
  AgeEffect_figure <- ggplot(axis, aes(x=SA.axis_rank, y=GAM.age.AdjRsq, fill = SA.axis_rank_signif)) + 
    geom_point(color = "gray", shape = 21, size=3.5) + 
    paletteer::scale_fill_paletteer_c("grDevices::RdYlBu", direction = -1, limits = c(min(axis$SA.axis_rank), max(axis$SA.axis_rank)), oob = squish) +
    paletteer::scale_color_paletteer_c("grDevices::RdYlBu", direction = -1, limits = c(min(axis$SA.axis_rank), max(axis$SA.axis_rank)), oob = squish) +
      labs(fill = "SA Axis Rank", x="\nSensorimotor-Association Axis Rank\n", y=expression(paste("Age Effect (Delta Adj", " R"^2, ")"))) +
      geom_smooth(data = axis, method='lm', se=TRUE, fill=alpha(c("gray70"),.9), col="black") +
      theme(
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.line = element_line(color = "black"),
        axis.text=element_text(size=24, color = "black"),
        panel.background=element_blank(),  
        legend.position = "bottom",
        legend.key.height = unit(1, 'cm'),
        legend.key.width = unit(2.3, 'cm'),
        legend.margin=margin(10,0,0,0),
        legend.text = element_text(size=24),
        legend.title = element_blank(),
        plot.title=element_text(hjust=0.5, size=24),
        plot.margin = margin(1, 1, 1, 1, "cm")) + scale_y_continuous(breaks = c(-0.05, 0, 0.05, 0.1, 0.15, 0.2), limits=c(ylim_lower, ylim_upper)) +     
     annotate(geom="text", x=x_pos, y=y_pos, label=r_text, color="black", size=7) + ggtitle(dataset)
  return(AgeEffect_figure)
}
 

 
r_text_GBC.schaefer200 <- expression(paste(italic("r"), " = -0.71, ",  , italic(p[spin]), "< 0.0001"))
GBC_schaefer200_PNC <- gam_figure3_outliers("PNC", GBC.axis.PNC, "GBC", "Schaefer200", r_text_GBC.schaefer200, 
                                      x_pos= 140, y_pos=0.15, ylim_lower=-0.09, ylim_upper=0.18)


r_text_GBC.schaefer200 <- expression(paste(italic("r"), " = -0.56, ",  , italic(p[spin]), "< 0.0001"))
GBC_schaefer200_NKI <- gam_figure3_outliers("NKI", GBC.axis.NKI, "GBC", "Schaefer200", r_text_GBC.schaefer200, 
                                      x_pos= 140, y_pos=0.15, ylim_lower=-0.09, ylim_upper=0.18)


r_text_GBC.schaefer200 <- expression(paste(italic("r"), " = -0.62, ",  , italic(p[spin]), "< 0.0001"))
GBC_schaefer200_HCPD <- gam_figure3_outliers("HCP-D", GBC.axis.HCPD, "GBC", "Schaefer200", r_text_GBC.schaefer200, 
                                  x_pos= 140, y_pos=0.15, ylim_lower=-0.09, ylim_upper=0.18)


r_text_GBC.schaefer200 <- expression(paste(italic("r"), " = -0.72 ",  , italic(p[spin]), "< 0.0001"))
GBC_schaefer200_HBN <- gam_figure3_outliers("HBN", GBC.axis.HBN, "GBC", "Schaefer200", r_text_GBC.schaefer200, 
                                  x_pos= 140, y_pos=0.15, ylim_lower=-0.09, ylim_upper=0.18)

 
figure3_alldataset <- ggarrange(GBC_schaefer200_PNC, GBC_schaefer200_HCPD, GBC_schaefer200_NKI, GBC_schaefer200_HBN, common.legend = TRUE, legend="bottom", nrow=2, ncol=2, labels=c("a", "c", "b", "d"), font.label = list(size = 24))

x.grob <- textGrob("S-A Axis Rank", 
                   gp=gpar(col="black", fontsize=28))
y.grob <- textGrob(expression(paste("Age Effect (Delta Adj", " R"^2, ")")), 
                   gp=gpar(col="black", fontsize=28), rot=90)

grid.arrange(arrangeGrob(figure3_alldataset, left = y.grob, bottom = x.grob))

fig3 <-arrangeGrob(figure3_alldataset, left = y.grob, bottom = x.grob)
ggsave("/Users/audluo/Library/CloudStorage/Box-Box/Box_PhD_Land/PennLINC/network_replication/Writing/Manuscript/FinalManuscript/FinalFigures/InDesign/Fig3.pdf", fig3, width=13, height=13, units="in")
write.csv(GBC.axis.PNC, "/cbica/projects/network_replication/manuscript/figure_data/figure3a_PNC_GBC_age_effect.csv")
write.csv(GBC.axis.NKI, "/cbica/projects/network_replication/manuscript/figure_data/figure3b_NKI_GBC_age_effect.csv")
write.csv(GBC.axis.HCPD, "/cbica/projects/network_replication/manuscript/figure_data/figure3c_HCPD_GBC_age_effect.csv")
write.csv(GBC.axis.HBN, "/cbica/projects/network_replication/manuscript/figure_data/figure3d_HBN_GBC_age_effect.csv")

Figure 4

gam.GBC.fitted.PNC <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/%2$s_%4$s.csv", derivs_root, "PNC", "GBC", "gam.GBC.fitted.schaefer200"))
gam.GBC.fitted.NKI <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/%2$s_%4$s.csv", derivs_root, "NKI", "GBC", "gam.GBC.fitted.schaefer200"))
gam.GBC.fitted.HCPD <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/%2$s_%4$s.csv", derivs_root, "HCPD", "GBC", "gam.GBC.fitted.schaefer200"))
gam.GBC.fitted.HBN <- read.csv(sprintf("%1$s%2$s/%3$s/GAM/%2$s_%4$s.csv", derivs_root, "HBN", "GBC", "gam.GBC.fitted.schaefer200"))

# load posterior fitted GBC correlations to SA-axis 
fits.SAaxis.posteriorcorrs.PNC <- readRDS(sprintf("%1$s%2$s/%3$s/GAM/%2$s_SAaxis_posteriorfits_correlation_byage_%4$s.RData", derivs_root, "PNC", "GBC", "schaefer200"))
colnames(fits.SAaxis.posteriorcorrs.PNC) <- sprintf("draw%s",seq(from = 1, to = ncol(fits.SAaxis.posteriorcorrs.PNC)))
fits.SAaxis.posteriorcorrs.PNC <- cbind(fits.SAaxis.posteriorcorrs.PNC, (gam.GBC.fitted.PNC %>% group_by(age) %>% group_keys()))
 
fits.SAaxis.posteriorcorrs.NKI <- readRDS(sprintf("%1$s%2$s/%3$s/GAM/%2$s_SAaxis_posteriorfits_correlation_byage_%4$s.RData", derivs_root, "NKI", "GBC", "schaefer200"))
colnames(fits.SAaxis.posteriorcorrs.NKI) <- sprintf("draw%s",seq(from = 1, to = ncol(fits.SAaxis.posteriorcorrs.NKI)))
fits.SAaxis.posteriorcorrs.NKI <- cbind(fits.SAaxis.posteriorcorrs.NKI, (gam.GBC.fitted.NKI %>% group_by(age) %>% group_keys()))
 
    
fits.SAaxis.posteriorcorrs.HCPD <- readRDS(sprintf("%1$s%2$s/%3$s/GAM/%2$s_SAaxis_posteriorfits_correlation_byage_%4$s.RData", derivs_root, "HCPD", "GBC", "schaefer200"))
colnames(fits.SAaxis.posteriorcorrs.HCPD) <- sprintf("draw%s",seq(from = 1, to = ncol(fits.SAaxis.posteriorcorrs.HCPD)))
fits.SAaxis.posteriorcorrs.HCPD <- cbind(fits.SAaxis.posteriorcorrs.HCPD, (gam.GBC.fitted.HCPD %>% group_by(age) %>% group_keys()))
 

fits.SAaxis.posteriorcorrs.HBN <- readRDS(sprintf("%1$s%2$s/%3$s/GAM/%2$s_SAaxis_posteriorfits_correlation_byage_%4$s.RData", derivs_root, "HBN", "GBC", "schaefer200"))
colnames(fits.SAaxis.posteriorcorrs.HBN) <- sprintf("draw%s",seq(from = 1, to = ncol(fits.SAaxis.posteriorcorrs.HBN)))
fits.SAaxis.posteriorcorrs.HBN <- cbind(fits.SAaxis.posteriorcorrs.HBN, (gam.GBC.fitted.HBN %>% group_by(age) %>% group_keys()))

compute_credible_intervals <- function(dataset) {
  fits.SAaxis.posteriorcorrs_GBC <- get(paste0("fits.SAaxis.posteriorcorrs.", dataset))
  gam.GBC.fitted.schaefer200 <- get(paste0("gam.GBC.fitted.", dataset))
  fits.SAaxis.posteriorcorrs <- fits.SAaxis.posteriorcorrs_GBC %>% select(contains("draw"))
  fits.SAaxis.mediancorr <- apply(fits.SAaxis.posteriorcorrs, 1, function(x){median(x)})
  cor.credible.intervals <- apply(fits.SAaxis.posteriorcorrs, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws 
  cor.credible.intervals <- as.data.frame(t(cor.credible.intervals))
    
  cor.credible.intervals <- cbind(cor.credible.intervals, (gam.GBC.fitted.schaefer200 %>% group_by(age) %>% group_keys()))
  cor.credible.intervals <- cbind(cor.credible.intervals, fits.SAaxis.mediancorr)
  colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")
      
  return(cor.credible.intervals)
}

cor.credible.intervals.PNC <- compute_credible_intervals("PNC")
cor.credible.intervals.NKI <- compute_credible_intervals("NKI")
cor.credible.intervals.HCPD <- compute_credible_intervals("HCPD")
cor.credible.intervals.HBN <- compute_credible_intervals("HBN")

plot_credible_intervals <- function(dataset) {
  cor.credible.intervals <- get(paste0("cor.credible.intervals.", dataset))
  ggplot(cor.credible.intervals, aes(x = age, y = abs(SA.correlation), ymin = abs(lower), ymax = abs(upper))) + 
    geom_line(size = .5) +
    geom_ribbon(alpha = .2, fill = c("grey60")) + geom_hline(yintercept=0.6,linetype=2) + 
    labs(x="Age", y="Spatial Alignment of GBC to S-A Axis (r)", title=dataset) + 
    theme_classic() + ylim(0, 0.62) + xlim(5, 23) + 
    theme(axis.text = element_text(size = 24, 
                                   color = c("black")), 
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          axis.line = element_line(size =.22), 
          axis.ticks = element_line(size = .22), 
          plot.title = element_text(size=24, hjust=0.5))  
 
}
 
PNC_ci_plot <- plot_credible_intervals("PNC") 
NKI_ci_plot <- plot_credible_intervals("NKI") 
HCPD_ci_plot <- plot_credible_intervals("HCPD") 
HBN_ci_plot <- plot_credible_intervals("HBN")
 
figure4_alldataset <- ggarrange(PNC_ci_plot, HCPD_ci_plot, NKI_ci_plot, HBN_ci_plot, labels='auto', font.label=list(size=20))

x.grob <- textGrob("Age (years)", 
                   gp=gpar(fontface="bold", col="black", fontsize=28))
y.grob <- textGrob("Absolute Spatial Alignment of FC Strength with S-A axis (r)", 
                   gp=gpar(fontface="bold", col="black", fontsize=28), rot=90)

grid.arrange(arrangeGrob(figure4_alldataset, left = y.grob, bottom = x.grob))

write.csv(cor.credible.intervals.PNC, "/cbica/projects/network_replication/manuscript/figure_data/figure4a_PNC_cor.credible.intervals_GBC_fitted.csv")
write.csv(cor.credible.intervals.NKI, "/cbica/projects/network_replication/manuscript/figure_data/figure4b_NKI_cor.credible.intervals_GBC_fitted.csv")
write.csv(cor.credible.intervals.HCPD, "/cbica/projects/network_replication/manuscript/figure_data/figure4c_HCPD_cor.credible.intervals_GBC_fitted.csv")
write.csv(cor.credible.intervals.HBN, "/cbica/projects/network_replication/manuscript/figure_data/figure4d_HBN_cor.credible.intervals_GBC_fitted.csv")

Figure 5

## FIGURE: SA axis average rank vs. age effect - with spin test
# included all datapoints
# @param axis A dataframe with SA-axis rank and GAM results
# @param metric A character string of connectivity metric (i.e. "GBC")
# @param atlas A character string of atlas name
# @param r_text An expression including info on p_spin and r
# @param x_pos A numeric for horizontal position of r_text
# @param y_pos A numeric for vertical position of r_text
gam_figure_p.spin <- function(axis, metric, atlas, r_text, x_pos, y_pos, ylim_lower, ylim_upper){
  axis <- axis %>% mutate(SA.axis_rank_signif = ifelse(significant.fdr == 1, SA.axis_rank, NA))
  AgeEffect_figure <- ggplot(axis, aes(x=SA.axis_rank, y=GAM.age.AdjRsq, fill = SA.axis_rank_signif)) + 
    geom_point(color = "gray", shape = 21, size=3.5) + 
    paletteer::scale_fill_paletteer_c("grDevices::RdYlBu", direction = -1, limits = c(min(axis$SA.axis_rank), max(axis$SA.axis_rank)), oob = squish) +
    paletteer::scale_color_paletteer_c("grDevices::RdYlBu", direction = -1, limits = c(min(axis$SA.axis_rank), max(axis$SA.axis_rank)), oob = squish) +
    ylim(ylim_lower, ylim_upper) +  # ylim(-0.065, 0.15) for GBC;  ylim(-0.07, 0.10) for BNC;  ylim(-0.055, 0.17) for WNC; 
      geom_smooth(data = axis, method='lm', se=TRUE, fill=alpha(c("gray70"),.9), col="black") +
      theme(
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.line = element_line(color = "black"),
        axis.text=element_text(size=24, color = "black"),
        panel.background=element_blank(),  
        legend.position = "none",
        plot.margin = margin(1, 1, 1, 1, "cm")) +
     annotate(geom="text", x=x_pos, y=y_pos, label=r_text, color="black", size=7)
  return(AgeEffect_figure)
}

# make df.dev.spin.BNC.<dataset> (formatted df's for applying spin test)
df.dev.spin.BNC.PNC <- prepDfSpin("schaefer200","PNC", "BNC")
df.dev.spin.BNC.NKI <- prepDfSpin("schaefer200","NKI", "BNC")
df.dev.spin.BNC.HCPD <- prepDfSpin("schaefer200","HCPD", "BNC")
df.dev.spin.BNC.HBN <- prepDfSpin("schaefer200","HBN", "BNC")


# make df.dev.spin.WNC.<dataset> (formatted df's for applying spin test)
df.dev.spin.WNC.PNC <- prepDfSpin("schaefer200","PNC", "WNC")
df.dev.spin.WNC.NKI <- prepDfSpin("schaefer200","NKI", "WNC")
df.dev.spin.WNC.HCPD <- prepDfSpin("schaefer200","HCPD", "WNC")
df.dev.spin.WNC.HBN <- prepDfSpin("schaefer200","HBN", "WNC")

# compute spatial correlation between age effect and S-A axis rank
PNC_corr_BNC <- round(cor.test(BNC.axis.PNC$GAM.age.AdjRsq, as.numeric(BNC.axis.PNC$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
NKI_corr_BNC <- round(cor.test(BNC.axis.NKI$GAM.age.AdjRsq, as.numeric(BNC.axis.NKI$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
HCPD_corr_BNC <- round(cor.test(BNC.axis.HCPD$GAM.age.AdjRsq, as.numeric(BNC.axis.HCPD$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
HBN_corr_BNC <- round(cor.test(BNC.axis.HBN$GAM.age.AdjRsq, as.numeric(BNC.axis.HBN$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)

PNC_corr_WNC <- round(cor.test(WNC.axis.PNC$GAM.age.AdjRsq, as.numeric(WNC.axis.PNC$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
NKI_corr_WNC <- round(cor.test(WNC.axis.NKI$GAM.age.AdjRsq, as.numeric(WNC.axis.NKI$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
HCPD_corr_WNC <- round(cor.test(WNC.axis.HCPD$GAM.age.AdjRsq, as.numeric(WNC.axis.HCPD$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)
HBN_corr_WNC <- round(cor.test(WNC.axis.HBN$GAM.age.AdjRsq, as.numeric(WNC.axis.HBN$SA.axis_rank), method=c("spearman"), exact=F)$estimate, 2)


perm.id.full_schaefer200 <- readRDS("/cbica/projects/network_replication/software/rotate_parcellation/schaefer200x7.coords_sphericalrotations_N10k_seed10.rds")
 
# spin permutation tests:
PNC_pspin_BNC <- perm.sphere.p(as.numeric(df.dev.spin.BNC.PNC$GAM.age.AdjRsq), as.numeric(df.dev.spin.BNC.PNC$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

NKI_pspin_BNC <- perm.sphere.p(as.numeric(df.dev.spin.BNC.NKI$GAM.age.AdjRsq), as.numeric(df.dev.spin.BNC.NKI$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

HCPD_pspin_BNC <- perm.sphere.p(as.numeric(df.dev.spin.BNC.HCPD$GAM.age.AdjRsq), as.numeric(df.dev.spin.BNC.HCPD$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

HBN_pspin_BNC <- perm.sphere.p(as.numeric(df.dev.spin.BNC.HBN$GAM.age.AdjRsq), as.numeric(df.dev.spin.BNC.HBN$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

corr_BNC_pspin <- data.frame(cbind(c(PNC_corr_BNC, NKI_corr_BNC, HCPD_corr_BNC, HBN_corr_BNC), c(PNC_pspin_BNC, NKI_pspin_BNC, HCPD_pspin_BNC, HBN_pspin_BNC)))  
corr_BNC_pspin <- cbind(c("PNC_BNC", "NKI_BNC", "HCPD_BNC", "HBN_BNC"), corr_BNC_pspin) %>% setNames(c("dataset_metric","rho", "pspin"))
rownames(corr_BNC_pspin) <- NULL

PNC_pspin_WNC <- perm.sphere.p(as.numeric(df.dev.spin.WNC.PNC$GAM.age.AdjRsq), as.numeric(df.dev.spin.WNC.PNC$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

NKI_pspin_WNC <- perm.sphere.p(as.numeric(df.dev.spin.WNC.NKI$GAM.age.AdjRsq), as.numeric(df.dev.spin.WNC.NKI$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')
  
HCPD_pspin_WNC <- perm.sphere.p(as.numeric(df.dev.spin.WNC.HCPD$GAM.age.AdjRsq), as.numeric(df.dev.spin.WNC.HCPD$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')

HBN_pspin_WNC <- perm.sphere.p(as.numeric(df.dev.spin.WNC.HBN$GAM.age.AdjRsq), as.numeric(df.dev.spin.WNC.HBN$SA.axis_rank),
                            perm.id.full_schaefer200, corr.type='spearman')  

corr_WNC_pspin <- data.frame(cbind(c(PNC_corr_WNC, NKI_corr_WNC, HCPD_corr_WNC, HBN_corr_WNC), c(PNC_pspin_WNC, NKI_pspin_WNC, HCPD_pspin_WNC, HBN_pspin_WNC)))  
corr_WNC_pspin <- cbind(c("PNC_WNC", "NKI_WNC", "HCPD_WNC", "HBN_WNC"), corr_WNC_pspin) %>% setNames(c("dataset_metric","rho", "pspin"))
rownames(corr_WNC_pspin) <- NULL
 
corr_pspin <- rbind(corr_BNC_pspin, corr_WNC_pspin)

write.csv(BNC.axis.PNC, "/cbica/projects/network_replication/manuscript/figure_data/figure5a_PNC_BNC.csv")
write.csv(WNC.axis.PNC, "/cbica/projects/network_replication/manuscript/figure_data/figure5e_PNC_WNC.csv")

write.csv(BNC.axis.NKI, "/cbica/projects/network_replication/manuscript/figure_data/figure5b_NKI_BNC.csv")
write.csv(WNC.axis.NKI, "/cbica/projects/network_replication/manuscript/figure_data/figure5f_NKI_WNC.csv")

write.csv(BNC.axis.HCPD, "/cbica/projects/network_replication/manuscript/figure_data/figure5c_HCPD_BNC.csv")
write.csv(WNC.axis.HCPD, "/cbica/projects/network_replication/manuscript/figure_data/figure5g_HCPD_WNC.csv")

write.csv(BNC.axis.HBN, "/cbica/projects/network_replication/manuscript/figure_data/figure5d_HBN_BNC.csv")
write.csv(WNC.axis.HBN, "/cbica/projects/network_replication/manuscript/figure_data/figure5h_HBN_WNC.csv")

r_text_BNC.PNC <- expression(paste(italic("r"), "= -0.66, ",  , italic(p[spin]), "< 0.0001"))
BNC_schaefer200_PNC <- gam_figure_p.spin(BNC.axis.PNC, "BNC", "schaefer200x7", r_text_BNC.PNC, 
                                        x_pos= 100, y_pos=0.1, ylim_lower=-0.07, ylim_upper=0.10)

r_text_WNC.PNC <- expression(paste(italic("r"), "= -0.49, ",  , italic(p[spin]), "< 0.001"))
WNC_schaefer200_PNC <- gam_figure_p.spin(WNC.axis.PNC, "WNC", "schaefer200x7", r_text_WNC.PNC, 
                                        x_pos= 100, y_pos=0.15, ylim_lower=-0.055, ylim_upper=0.17)

r_text_BNC.NKI <- expression(paste(italic("r"), "= -0.50, ",  , italic(p[spin]), "< 0.001"))
BNC_schaefer200_NKI <- gam_figure_p.spin(BNC.axis.NKI, "BNC", "schaefer200x7", r_text_BNC.NKI, 
                                        x_pos= 100, y_pos=0.1, ylim_lower=-0.07, ylim_upper=0.10)

r_text_WNC.NKI <- expression(paste(italic("r"), "= -0.31, ",  , italic(p[spin]), "< 0.05"))
WNC_schaefer200_NKI <- gam_figure_p.spin(WNC.axis.NKI, "WNC", "schaefer200x7", r_text_WNC.NKI, 
                                        x_pos= 100, y_pos=0.16, ylim_lower=-0.055, ylim_upper=0.17)


r_text_BNC.HCPD <- expression(paste(italic("r"), "= -0.55, ",  , italic(p[spin]), "< 0.0001"))
BNC_schaefer200_HCPD <- gam_figure_p.spin(BNC.axis.HCPD, "BNC", "schaefer200x7", r_text_BNC.HCPD, 
                                        x_pos= 100, y_pos=0.1, ylim_lower=-0.07, ylim_upper=0.10)

r_text_WNC.HCPD <- expression(paste(italic("r"), "= -0.31, ",  , italic(p[spin]), "< 0.05"))
WNC_schaefer200_HCPD <- gam_figure_p.spin(WNC.axis.HCPD, "WNC", "schaefer200x7", r_text_WNC.HCPD, 
                                        x_pos= 100, y_pos=0.16, ylim_lower=-0.055, ylim_upper=0.17)


r_text_BNC.HBN <- expression(paste(italic("r"), "= -0.70, ",  , italic(p[spin]), "< 0.0001"))
BNC_schaefer200_HBN <- gam_figure_p.spin(BNC.axis.HBN, "BNC", "schaefer200x7", r_text_BNC.HBN, 
                                        x_pos= 100, y_pos=0.1, ylim_lower=-0.07, ylim_upper=0.10)

r_text_WNC.HBN <- expression(paste(italic("r"), "= -0.38, ",  , italic(p[spin]), "< 0.001"))
WNC_schaefer200_HBN <- gam_figure_p.spin(WNC.axis.HBN, "WNC", "schaefer200x7", r_text_WNC.HBN, 
                                        x_pos= 100, y_pos=0.15, ylim_lower=-0.055, ylim_upper=0.17)


figure5_alldatasets <- ggarrange(BNC_schaefer200_PNC + labs(title = "Between-Network Connectivity") + 
                                 theme(plot.title=element_text(size=20, hjust=0.5, vjust=5)), 
          WNC_schaefer200_PNC + labs(title = "Within-Network Connectivity") + 
                                 theme(plot.title=element_text(size=20, hjust=0.5, vjust=5)), 
          BNC_schaefer200_NKI, WNC_schaefer200_NKI, 
          BNC_schaefer200_HCPD, WNC_schaefer200_HCPD, 
          BNC_schaefer200_HBN, WNC_schaefer200_HBN, ncol=2, nrow=4, labels=c("a", "e", "b", "f", "c", "g", "d", "h"))

x.grob <- textGrob("S-A Axis Rank", 
                   gp=gpar(col="black", fontsize=28))
y.grob <- textGrob(expression(paste("Age Effect (Delta Adj", " R"^2, ")")), 
                   gp=gpar(col="black", fontsize=28), rot=90)

grid.arrange(arrangeGrob(figure5_alldatasets, left = y.grob, bottom = x.grob))

# row 1 = PNC
# row 2 = NKI
# row 3 = HCP-D
# row 4 = HBN

Figure 6

# model the surface
make_surface_plot_datasets <- function(dataset) {
  double_age_SA.diff <- get(paste0("double_age_SA.diff_", dataset))
  g2 <- gam(GAM.age.AdjRsq ~ te(parcel2.SA.vec,parcel1.SA.vec,k=3),data = double_age_SA.diff)
 
  surface.plot.fig <- gg_tensor(g2)+theme_classic() +
    theme(
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      axis.line = element_line(color = "black"),
      axis.text=element_text(size=24, color = "black")) + 
   scale_fill_gradientn(name = expression(paste("Age Effect (Delta Adj", " R"^2, ")")), colors= c("#8a0f63", "#FFFFFF","#f4b100"), 
                        limits = c(-0.014,0.015), oob=squish, values=rescale(c(-0.014,0,0.015))) +
    geom_vline(xintercept = mean(range(double_age_SA.diff$parcel1.SA.vec)),
               linetype='dashed',size=1) + 
    geom_hline(yintercept = mean(range(double_age_SA.diff$parcel1.SA.vec)),
               linetype='dashed',size=1) + labs(title=dataset) + 
    theme(legend.position='bottom',
          legend.key.height = unit(1, 'cm'),
          legend.key.width = unit(3, 'cm'),
          legend.title = element_text(size=20),
          legend.text = element_text(size=20),
          strip.background = element_blank(),
          plot.margin = margin(1, 1, 1, 1, "cm"),
          plot.title = element_text(size=28, hjust=0.5))  
  return(surface.plot.fig)
}
 
 
surfacePlot_PNC <- make_surface_plot_datasets("PNC")
surfacePlot_NKI <- make_surface_plot_datasets("NKI")
surfacePlot_HCPD <- make_surface_plot_datasets("HCPD")
surfacePlot_HBN <- make_surface_plot_datasets("HBN")

figure6_alldatasets<- ggarrange(surfacePlot_PNC, surfacePlot_HCPD, surfacePlot_NKI, surfacePlot_HBN, common.legend = TRUE, labels = c("a", "c", "b", "d"))
 

x.grob <- textGrob("S-A Axis Rank", 
                   gp=gpar(fontface="bold", col="black", fontsize=28))
y.grob <- textGrob("S-A Axis Rank", 
                   gp=gpar(fontface="bold", col="black", fontsize=28), rot=90)


grid.arrange(arrangeGrob(figure6_alldatasets, left = y.grob, bottom = x.grob))

write.csv(double_age_SA.diff_PNC, "/cbica/projects/network_replication/manuscript/figure_data/figure6a_PNC_edge.csv")
write.csv(double_age_SA.diff_NKI, "/cbica/projects/network_replication/manuscript/figure_data/figure6b_NKI_edge.csv")
write.csv(double_age_SA.diff_HCPD, "/cbica/projects/network_replication/manuscript/figure_data/figure6c_HCPD_edge.csv")
write.csv(double_age_SA.diff_HBN, "/cbica/projects/network_replication/manuscript/figure_data/figure6d_HBN_edge.csv")